
#ifndef _TRILINEAR_H
#define _TRILINEAR_H_

#include "Tuna.hpp"

class Trilinear 
{
///
///  Traits to define the type of arrays to be used.
///
///  floatTinyArray_t : Array of floats mostly used for deltas of meshes
///    
  typedef TunaArray<double, 3>::tiny floatTinyArray_t;
  typedef TunaArray<double, 3>::huge ScalarField3D;
 
 public:
    Trilinear() { };
    ~Trilinear() { };

/*!
 * This function locate the position of the particle on the 3D mesh. 
 * The cell where the particle is located is identified with the indexes 
 * (I, J, K). The local position inside the cell (I,J,K) is defined by 
 * (alpha, beta, gamma).
 * \param pos Position of the particle.
 * \param deltas Deltas on each axis of the mesh.
 */
    inline void locate(floatTinyArray_t &pos, const floatTinyArray_t &deltas)
	{
	    double a = pos(0) / deltas (0);
	    double b = pos(1) / deltas (1);
	    double c = pos(2) / deltas (2);
	    
	    I = static_cast<int>(a);
	    J = static_cast<int>(b);
	    K = static_cast<int>(c);
	    
	    alpha = a - I;
	    beta  = b - J;
	    gamma = c - K;
	}
	
/*!
 * This is a trilinear interpolation of the velocity to the actual position
 * of the particle. The result is stored in the array UVW.
 * \param pos Position of the particle.
 * \param deltas Deltas on each axis of the mesh.
 * \param u X-component of the velocity.
 * \param v Y-component of the velocity.
 * \param w Z-component of the velocity.
 */  
    inline floatTinyArray_t interpolate(floatTinyArray_t &pos,
					const floatTinyArray_t &deltas,
					ScalarField3D &u,
					ScalarField3D &v,
					ScalarField3D &w)
	{
	    locate(pos, deltas);
	    UVW(0) = ( ( u(I  ,J  ,K  ) * (1 - alpha) +
			 u(I+1,J  ,K  ) * alpha ) * (1 - beta) +
		       ( u(I  ,J+1,K  ) * (1 - alpha) +
			 u(I+1,J+1,K  ) * alpha ) * beta ) * (1 - gamma) +
		     ( ( u(I  ,J  ,K+1) * (1 - alpha) +
			 u(I+1,J  ,K+1) * alpha ) * (1 - beta) +
		       ( u(I  ,J+1,K+1) * (1 - alpha ) +
			 u(I+1,J+1,K+1) * alpha ) * beta ) * gamma;
	    
	    UVW(1) = ( ( v(I  ,J  ,K  ) * (1 - alpha) +
			 v(I+1,J  ,K  ) * alpha ) * (1 - beta) +
		       ( v(I  ,J+1,K  ) * (1 - alpha) +
			 v(I+1,J+1,K  ) * alpha ) * beta ) * (1 - gamma) +
		     ( ( v(I  ,J  ,K+1) * (1 - alpha) +
			 v(I+1,J  ,K+1) * alpha ) * (1 - beta) +
		       ( v(I  ,J+1,K+1) * (1 - alpha ) +
			 v(I+1,J+1,K+1) * alpha ) * beta ) * gamma;
		     
	    UVW(2) = ( ( w(I  ,J  ,K  ) * (1 - alpha) +
			 w(I+1,J  ,K  ) * alpha ) * (1 - beta) +
		       ( w(I  ,J+1,K  ) * (1 - alpha) +
			 w(I+1,J+1,K  ) * alpha ) * beta ) * (1 - gamma) +
		     ( ( w(I  ,J  ,K+1) * (1 - alpha) +
			 w(I+1,J  ,K+1) * alpha ) * (1 - beta) +
		       ( w(I  ,J+1,K+1) * (1 - alpha ) +
			 w(I+1,J+1,K+1) * alpha ) * beta ) * gamma;
	    
	    return UVW;
	}

 private:
    int I, J, K ;
    double alpha, beta, gamma;  
    floatTinyArray_t UVW;
};

#endif // _TRILINEAR_H_
